function res = showSensitivity(obj, channelIndex, sampleName, calSigFreq, calSigAmpl, varargin)
% ZXD version of analyze() method.
% channelIndex, sampleName : see getNodeData()
% calSigFreq, calSigAmpl: the frequency and amplitude of the small magnetic field signal for calibration. In unit of Hz and nTpk. The unit 'nT' can be changed by parameter "unit_calSigAmpl"
    spect = obj.getNodeData(channelIndex, sampleName);
    freq = spect.freq(:).';
    val = spect.value(:).';
    header = spect.info;
    df = freq(2) - freq(1);
    
    p=inputParser;
    p.addParameter('func', @semilogy, @(x) ismember(x, {@semilogy, @plot, @loglog}));
    p.addParameter('uBound', max(val), @(x) isnumeric(x) && isscalar(x) );             % upperbound of histcounts, can be used to exclude peaks
    p.addParameter('nBin', 100, @(x) isscalar(x) && x > 1);                            % nBin for histcounts
    p.addParameter('unit_calSigAmpl', 'nT', @(x) ischar(x));                           % unit of magnetic field
    p.addParameter('search_range', max(10*df, 1), @(x) isscalar(x) && x>=df);          % the freq. range to search for calibration signal peak.
    p.addParameter('range', minmax(freq), @(x) length(x) == 2 && x(2) > x(1) );        % freq. range for noise background statistics 
    p.addParameter('isPlot', true, @islogical);
    p.addParameter('figHdl', [], @(x) isempty(x) || isa(x, 'matlab.ui.Figure'));
    p.parse(varargin{:});
    
    magUnit = p.Results.unit_calSigAmpl;
    
    %%
    enbw = header.nenbw*df;

    rangeIdx = ( freq >= calSigFreq - p.Results.search_range & freq <= calSigFreq + p.Results.search_range );
    [pks, res.location] = findpeaks(val(rangeIdx), freq(rangeIdx), 'Npeak', 1, 'SortStr', 'descend');
    res.Vpeak = pks * sqrt(enbw);
    % pks:       spectral density, in unit of Vpk / sqrt(Hz)
    % res.Vpeak: FFT spectrum, in unit of Vpk
    
    rangeIdx = ( freq >= p.Results.range(1) & freq <= p.Results.range(2) );
	uBound = p.Results.uBound * 2;
	res.stdBg = 0;
	binSize = Inf;
	while res.stdBg<binSize*p.Results.nBin/20
		uBound = uBound/2;
		[c, edges]=histcounts(val(val<uBound & rangeIdx), p.Results.nBin);
		binSize = edges(2) - edges(1);
		edgesVal = edges(1:end-1)+0.5*binSize;
		[res.background, res.stdBg] = fitGaussian(edgesVal, c);
	end
    
    % conversion factor: convert voltage to magnetic field
    factor = calSigAmpl/res.Vpeak;
    res.sensitivity = factor * res.background;
    res.stdSen = factor * res.stdBg;
    res.unitSen = [magUnit '/sqrt(Hz)'];
    res.gain = 1e-3*factor;
    res.unitGain = [magUnit '/mV'];

    %%
    if p.Results.isPlot
        if isempty(p.Results.figHdl)
            figure('Name', 'Spectrum Analyser', 'Position', [100, 100, 1200,800]);
        else
            gcf = p.Results.figHdl;
        end
        ax1 = subplot(1, 4, 1); ax2 = subplot(1, 4, [2 3 4]);
        semilogx(ax1, edgesVal*factor, c, '-', [res.sensitivity, res.sensitivity], minmax(c), 'r--'); grid(ax1, 'on');
        p.Results.func(ax2, freq, val*factor, minmax(freq), [res.sensitivity, res.sensitivity], 'r--'); grid(ax2, 'on'); hold(ax2, 'on');
        
        if p.Results.range(1) > min(freq)
            p.Results.func(ax2, [p.Results.range(1) p.Results.range(1)], ylim(ax2), 'k--');
        end
        if p.Results.range(2) < max(freq)            
            p.Results.func(ax2, [p.Results.range(2) p.Results.range(2)], ylim(ax2), 'k-.');
        end
        
        view(ax1, -90, 90);
        xlim(ax1, ylim(ax2)); xlabel(ax1, ['PSD (' magUnit '\cdot Hz^{-1/2})']); ylabel(ax1, 'Counts');
        
        
        xlim(ax2, minmax(freq));
        yticklabels(ax2, {});
        xlabel(ax2, 'Freuquency (Hz)');
        
        for k = 1:length(res.Vpeak)
            text(ax2, res.location(k), res.Vpeak(k)/sqrt(enbw)*factor, ...
                sprintf(['\\leftarrow Peak #%d = %.1f ' magUnit 'pk @ %.1fHz'], k, calSigAmpl, res.location),...
                'Color','red','FontSize',16);
        end
        text(ax2, mean(minmax(freq)), res.sensitivity+10.0*res.stdSen, ...
                  [sprintf('(%3.2e\\pm%2.1e) ',res.sensitivity, res.stdSen) magUnit ' \cdot Hz^{-1/2})'],...
                            'Color','red','FontSize',16);
        title(ax2,['Gain: ',num2str(res.gain,3),res.unitGain]);
    end
end


function [mu, sigma, amplitude, gof] = fitGaussian(x, y)
    [maxY, maxIdx] = max(y);
    x0 = x(maxIdx);
    [~, idx] = min(abs(y-maxY/exp(1.0)));
    stdGuess = abs(x(max(idx)) - x0);
    
    [xData, yData] = prepareCurveData( x, y );

    % Set up fittype and options.
    ft = fittype( 'gauss1' );  % a1*exp(-((x-b1)/c1)^2)
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Lower =      [0.9*maxY 0.5*x0 0];
    opts.Upper =      [1.1*maxY 2.0*x0 Inf];
    opts.StartPoint = [    maxY     x0 stdGuess];

    % Fit model to data.
    [fitresult, gof] = fit( xData, yData, ft, opts );
    mu = fitresult.b1;
    sigma = fitresult.c1/sqrt(2.0);
    amplitude = fitresult.a1;
end